#######################################################
#######################################################
#######################################################
### CREATED BY JONATHAN KING AND JESSICA SCHOENHERR ###
##### REPLICATION DATA FOR "A MATTER OF OPINION?" #####
#######################################################
#######################################################
#######################################################


library(readxl)
library(tidyverse)
library(ggeffects)
library(stargazer)
library(poliscidata)
library(ggplot2)
library(psych)
library(broom)
library(emmeans)

dataDeathPenalty <- read.csv("DeathPenaltyDataPart2-20221022.csv")

####################
### DATA SORTING ###
####################

dataDeathPenalty$treatmentGroup <- as.factor(dataDeathPenalty$treatmentGroup)
dataDeathPenalty$decisionTerm <- as.numeric(dataDeathPenalty$decisionTherm)
dataDeathPenalty$partisanship <- as.factor(dataDeathPenalty$partisanship)
dataDeathPenalty$female <- as.factor(dataDeathPenalty$female)

dataDeathPenalty <- within(dataDeathPenalty, treatmentGroup <- relevel(treatmentGroup, ref = 5))

dataDeathPenalty$partisanship <- sjmisc::rec(dataDeathPenalty$partisanship, rec = "1=Democrat;2=Independent;3=Republican", as.num = FALSE)

dataDeathPenalty$scTherm <- as.numeric(dataDeathPenalty$scTherm)
dataDeathPenalty$decisionTherm <- as.numeric(dataDeathPenalty$decisionTherm)

####################
### DESCRIPTIVES ###
####################

dataDeathPenalty %>% group_by(treatmentGroup, partisanship, female) %>% summarize(mean = mean(decisionTherm)) %>% print(n = 30)

dataDeathPenalty$count <- 1

dataDeathPenalty %>% group_by(treatmentGroup, partisanship, female) %>% count(count) %>% print(n = 30)

dataDeathPenalty %>% group_by(partisanship) %>% summarize(mean = mean(decisionTherm))
summary(aov(decisionTherm ~ partisanship, data = dataDeathPenalty))
TukeyHSD(aov(decisionTherm ~ partisanship, data = dataDeathPenalty))
dpPart <- lm(decisionTherm ~ partisanship, data = dataDeathPenalty)
car::Anova(dpPart)

dataDeathPenalty %>% group_by(female) %>% summarize(mean = mean(decisionTherm))
summary(aov(decisionTherm ~ female, data = dataDeathPenalty))
TukeyHSD(aov(decisionTherm ~ female, data = dataDeathPenalty))
dpFem <- lm(decisionTherm ~ female, data = dataDeathPenalty)
car::Anova(dpFem)

dataDeathPenalty %>% group_by(treatmentGroup) %>% summarize(mean = mean(decisionTherm))

################
### BASELINE ###
################

deathPenaltyTreatment <- lm(decisionTherm ~ treatmentGroup, 
						data = dataDeathPenalty)
summary(deathPenaltyTreatment)
nobs(deathPenaltyTreatment)

stargazer(deathPenaltyTreatment, type = "text")

car::Anova(deathPenaltyTreatment)

#dpCoef <- tidy(deathPenaltyTreatment)
#dpCoef$upperCI <- (dpCoef$std.error * 1.96) + dpCoef$estimate
#dpCoef$lowerCI <- (dpCoef$std.error * -1.96) + dpCoef$estimate

dpTreatmentPreds <- emmeans::emmeans(deathPenaltyTreatment, specs = "treatmentGroup")
pairs(dpTreatmentPreds)
dpTreatmentPreds <- as.data.frame(dpTreatmentPreds)

dpTreatmentPlot <- ggplot(dpTreatmentPreds) +
	geom_bar( aes(x = treatmentGroup, y = emmean), stat="identity", fill = "grey75") +
	geom_errorbar( aes(x = treatmentGroup, ymin = lower.CL, ymax = upper.CL), width=0.1) +
	theme_bw() +
	xlab("\nOpinion Writer Treatment Group")  +
	ylab("Participants' Warmth Toward Decision\n") +
	ggtitle("Death Penalty Vignette - Direct Treatment Effects") +
	theme(plot.title = element_text(hjust = 0.5)) +
	scale_x_discrete(labels = c("Control", "Liberal\nMale\nJustice", "Conservative\nMale\nJustice", "Liberal\nFemale\nJustice", "Conservative\nFemale\nJustice")) +
	scale_y_continuous(limits = c(0, 80), breaks = seq(0, 80, 10))
dpTreatmentPlot

#dpTreatmentCoef <- ggplot(dpCoef[2:5,], aes(y = estimate, x = term, ymin = lowerCI, ymax = upperCI)) + 
#	geom_pointrange(fatten = 2, size = 1) + 
#	theme_bw() +
#	xlab("\nOpinion Writer Treatment Group") +
#	ylab("Difference in Support Compared to Control\n") +
#	ggtitle("Death Penalty Vignette Coefficient Plot - Baseline Model") +
#	theme(plot.title = element_text(hjust = 0.5)) +
#	scale_x_discrete(labels = c("Liberal\nMale\nJustice", "Conservative\nMale\nJustice", "Liberal\nFemale\nJustice", "Conservative\nFemale\nJustice")) + 
#	scale_y_continuous(limits = c(-20, 10), breaks = seq(-20, 10, 5)) +
#	theme(legend.position = "none") +
#	geom_hline(yintercept = 0, size = 1, linetype = "dashed")
#dpTreatmentCoef

################################
### CONTROL FOR PARTISANSHIP ###
################################

deathPenaltyPartisanship <- lm(decisionTherm ~ treatmentGroup * partisanship,
						data = dataDeathPenalty)
summary(deathPenaltyPartisanship)
nobs(deathPenaltyPartisanship)

car::Anova(deathPenaltyPartisanship)

deathPenaltyPartisanshipTestDem <- data.frame(treatmentGroup = c(factor(6), factor(7), factor(8), factor(9), factor(10)), partisanship = "Democrat")
predict(deathPenaltyPartisanship, newdata = deathPenaltyPartisanshipTestDem)

deathPenaltyPartisanshipTestInd <- data.frame(treatmentGroup = c(factor(6), factor(7), factor(8), factor(9), factor(10)), partisanship = "Independent")
predict(deathPenaltyPartisanship, newdata = deathPenaltyPartisanshipTestInd)

deathPenaltyPartisanshipTestRep <- data.frame(treatmentGroup = c(factor(6), factor(7), factor(8), factor(9), factor(10)), partisanship = "Republican")
predict(deathPenaltyPartisanship, newdata = deathPenaltyPartisanshipTestRep)

deathPenaltyPartisanshipPreds <- emmeans::emmeans(deathPenaltyPartisanship, specs = "partisanship", by = "treatmentGroup")
pairs(deathPenaltyPartisanshipPreds)

deathPenaltyPartisanshipPredsDF <- as.data.frame(deathPenaltyPartisanshipPreds)

deathPenaltyPartisanPlot <- ggplot(deathPenaltyPartisanshipPredsDF, aes(y = emmean, x = treatmentGroup, color = partisanship, shape = partisanship)) +
	geom_point(position = position_dodge(0.3)) + 
	geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0, position = position_dodge(0.3)) +
	theme_bw() +
	xlab("\nOpinion Writer Treatment Group") +
	ylab("Predicted Value - Decision Thermometer\n") +
	ggtitle("Death Penalty Vignettes by Participant Partisanship") +
	theme(plot.title = element_text(hjust = 0.5)) +
	scale_x_discrete(labels = c("Control", "Liberal\nMale\nJustice", "Conservative\nMale\nJustice", "Liberal\nFemale\nJustice", "Conservative\nFemale\nJustice")) + 
	scale_y_continuous(limits = c(20, 100), breaks = seq(20, 100, 10)) +
	theme(legend.position = "bottom", legend.title = element_blank()) +
	scale_color_manual(values = c("black", "grey65", "grey50"), labels = c("Democrat", "Independent", "Republican")) + 
	scale_shape_manual(values = c(15, 17, 19), labels = c("Democrat", "Independent", "Republican"))
deathPenaltyPartisanPlot

deathPenaltyPartisanshipPreds2 <- emmeans::emmeans(deathPenaltyPartisanship, specs = "treatmentGroup", by = "partisanship")
pairs(deathPenaltyPartisanshipPreds2)

##########################
### CONTROL FOR GENDER ###
##########################

deathPenaltyGender <- lm(decisionTherm ~ treatmentGroup * female, 
				data = dataDeathPenalty)
summary(deathPenaltyGender)
nobs(deathPenaltyGender)

car::Anova(deathPenaltyGender)

deathPenaltyGenderTestMale <- data.frame(treatmentGroup = c(factor(6), factor(7), factor(8), factor(9), factor(10)), female = factor(0))
predict(deathPenaltyGender, newdata = deathPenaltyGenderTestMale)

deathPenaltyGenderTestFemale <- data.frame(treatmentGroup = c(factor(6), factor(7), factor(8), factor(9), factor(10)), female = factor(1))
predict(deathPenaltyGender, newdata = deathPenaltyGenderTestFemale)

deathPenaltyGenderPreds <- emmeans::emmeans(deathPenaltyGender, specs = "female", by = "treatmentGroup")
pairs(deathPenaltyGenderPreds)

deathPenaltyGenderPredsDF <- as.data.frame(deathPenaltyGenderPreds)

deathPenaltyGenderPreds2 <- emmeans::emmeans(deathPenaltyGender, specs = "treatmentGroup", by = "female")
pairs(deathPenaltyGenderPreds2)

deathPenaltyGenderPlot <- ggplot(deathPenaltyGenderPredsDF, aes(y = emmean, x = treatmentGroup, color = female, shape = female)) +
	geom_point(position = position_dodge(0.3)) + 
	geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0, position = position_dodge(0.3)) +
	theme_bw() +
	xlab("\nOpinion Writer Treatment Group") +
	ylab("Predicted Value - Decision Thermometer\n") +
	ggtitle("Death Penalty Vignettes by Participant Gender") +
	theme(plot.title = element_text(hjust = 0.5)) +
	scale_x_discrete(labels = c("Control", "Liberal\nMale\nJustice", "Conservative\nMale\nJustice", "Liberal\nFemale\nJustice", "Conservative\nFemale\nJustice")) + 
	scale_y_continuous(limits = c(20, 100), breaks = seq(20, 100, 10)) +
	theme(legend.position = "bottom", legend.title = element_blank()) +
	scale_color_manual(values = c("black", "grey50"), labels = c("Male", "Female")) +
	scale_shape_manual(values = c(15, 19), labels = c("Male", "Female"))
deathPenaltyGenderPlot

###############################
### GENDER AND PARTISANSHIP ###
###############################

modelDeathPenalty <- lm(decisionTherm ~ treatmentGroup
						+ female
						+ partisanship
						+ female * partisanship
						+ treatmentGroup * female * partisanship
						+ treatmentGroup * female
						+ treatmentGroup * partisanship, 
						data = dataDeathPenalty)
summary(modelDeathPenalty)
nobs(modelDeathPenalty)

stargazer(modelDeathPenalty, type = "text")

car::Anova(modelDeathPenalty)

deathPenaltyMaleDem <- data.frame(treatmentGroup = c(factor(6), factor(7), factor(8), factor(9), factor(10)), female = factor(0), partisanship = "Democrat")
predict(modelDeathPenalty, newdata = deathPenaltyMaleDem)

deathPenaltyMaleInd <- data.frame(treatmentGroup = c(factor(6), factor(7), factor(8), factor(9), factor(10)), female = factor(0), partisanship = "Independent")
predict(modelDeathPenalty, newdata = deathPenaltyMaleInd)

deathPenaltyMaleRep <- data.frame(treatmentGroup = c(factor(6), factor(7), factor(8), factor(9), factor(10)), female = factor(0), partisanship = "Republican")
predict(modelDeathPenalty, newdata = deathPenaltyMaleRep)

deathPenaltyFemaleDem <- data.frame(treatmentGroup = c(factor(6), factor(7), factor(8), factor(9), factor(10)), female = factor(1), partisanship = "Democrat")
predict(modelDeathPenalty, newdata = deathPenaltyFemaleDem)

deathPenaltyFemaleInd <- data.frame(treatmentGroup = c(factor(6), factor(7), factor(8), factor(9), factor(10)), female = factor(1), partisanship = "Independent")
predict(modelDeathPenalty, newdata = deathPenaltyFemaleInd)

deathPenaltyFemaleRep <- data.frame(treatmentGroup = c(factor(6), factor(7), factor(8), factor(9), factor(10)), female = factor(1), partisanship = "Republican")
predict(modelDeathPenalty, newdata = deathPenaltyFemaleRep)

deathPenaltyPreds <- emmeans::emmeans(modelDeathPenalty, specs = c("female", "partisanship"), by = "treatmentGroup")
pairs(deathPenaltyPreds)

deathPenaltyPreds <- emmeans::emmeans(modelDeathPenalty, specs = "treatmentGroup", by = c("female", "partisanship"))
pairs(deathPenaltyPreds)

deathPenaltyPredsDF <- as.data.frame(deathPenaltyPreds)

deathPenaltyPredsDF$group <- ifelse(deathPenaltyPredsDF$partisanship == "Democrat" & deathPenaltyPredsDF$female == 0, 1, 
ifelse(deathPenaltyPredsDF$partisanship == "Democrat" & deathPenaltyPredsDF$female == 1, 2,
ifelse(deathPenaltyPredsDF$partisanship == "Independent" & deathPenaltyPredsDF$female == 0, 3,
ifelse(deathPenaltyPredsDF$partisanship == "Independent" & deathPenaltyPredsDF$female == 1, 4,
ifelse(deathPenaltyPredsDF$partisanship == "Republican" & deathPenaltyPredsDF$female == 0, 5,
ifelse(deathPenaltyPredsDF$partisanship == "Republican" & deathPenaltyPredsDF$female == 1, 6, 0))))))
deathPenaltyPredsDF$group <- as.factor(deathPenaltyPredsDF$group)

#library(jtools)
#plots_coefs(modelDeathPenalty)

#dpNames <- as_labeller(c(`1`="Male Democrat Participants", `2`="Female Democrat Participants", `3`="Male Independent Participants", `4`="Female Independent Participants", `5`="Male Republican Participants", `6`="Female Republican Participants"))

#deathPenaltyPlot1 <- ggplot(deathPenaltyPredsDF, aes(y = emmean, x = treatmentGroup, ymin = lower.CL, max = upper.CL)) +
#	geom_pointrange(fatten = 2, size = 1) + 
#	geom_point(position = position_dodge(0.3)) + 
#	facet_wrap(~ group, scales = "free", labeller = dpNames) + 
#	theme_bw() +
#	xlab("\nOpinion Writer Treatment Group") +
#	ylab("Predicted Value - Decision Thermometer\n") +
#	ggtitle("Death Penalty Vignette Full Results") +
#	theme(plot.title = element_text(hjust = 0.5)) +
#	scale_x_discrete(labels = c("Control", "Liberal\nMale\nJustice", "Conservative\nMale\nJustice", "Liberal\nFemale\nJustice", "Conservative\nFemale\nJustice")) + 
#	scale_y_continuous(limits = c(20, 100), breaks = seq(20, 100, 20)) +
#	theme(legend.position = "bottom", legend.title = element_blank())
#deathPenaltyPlot1

#deathPenaltyPreds2 <- emmeans::emmeans(modelDeathPenalty, specs = c("treatmentGroup", "female"), by = "partisanship")
#pairs(deathPenaltyPreds2)

#deathPenaltyPredsDF$treatmentGroup <- factor(deathPenaltyPredsDF$treatmentGroup, levels = c(10, 6, 7, 8, 9))

#dpNames2 <- as_labeller(c(`10`="Control", `6`="Liberal Male Justice", `7`="Conservative Male Justice", `8`="Liberal Female Justice", `9`="Conservative Female Justice"))

#deathPenaltyPlot2 <- ggplot(deathPenaltyPredsDF, aes(y = emmean, x = partisanship, ymin = lower.CL, ymax = upper.CL, color = female, shape = female)) +
#	geom_pointrange(fatten = 2, size = 1, position = position_dodge(0.4)) + 
#	facet_grid(~treatmentGroup, labeller = dpNames2) + 
#	theme_bw() +
#	xlab("\nParticipant Partisanship") +
#	ylab("Predicted Value - Decision Thermometer\n") +
#	ggtitle("Death Penalty Vignette Expanded Results") +
#	theme(plot.title = element_text(hjust = 0.5)) +
#	scale_x_discrete(labels = c("Dem", "Ind", "Rep")) + 
#	scale_y_continuous(limits = c(20, 100), breaks = seq(20, 100, 10)) +
#	theme(legend.position = "bottom", legend.title = element_blank()) +
#	scale_color_manual(values = c("black", "grey50"), labels = c("Male", "Female")) +
#	scale_shape_manual(values = c(15, 19), labels = c("Male", "Female"))
#deathPenaltyPlot2

dpNames3 <- as_labeller(c(`Democrat` = "Democrat Participants", `Independent` = "Independent Participants", `Republican` = "Republican Participants"))

#deathPenaltyPlot3 <- ggplot(deathPenaltyPredsDF, aes(y = emmean, x = treatmentGroup, ymin = lower.CL, ymax = upper.CL, color = female, shape = female)) +
#	geom_pointrange(fatten = 2, size = 1, position = position_dodge(0.4)) + 
#	facet_grid(~partisanship, labeller = dpNames3) + 
#	theme_bw() +
#	xlab("\nOpinion Writer Treatment Groups") +
#	ylab("Predicted Value - Decision Thermometer\n") +
#	ggtitle("Death Penalty Vignette Expanded Results") +
#	theme(plot.title = element_text(hjust = 0.5)) +
#	scale_x_discrete(labels = c("Control", "Liberal\nMale\nJustice", "Conservative\nMale\nJustice", "Liberal\nFemale\nJustice", "Conservative\nFemale\nJustice")) + 
#	scale_y_continuous(limits = c(10, 100), breaks = seq(10, 100, 10)) +
#	theme(legend.position = "bottom", legend.title = element_blank()) +
#	scale_color_manual(values = c("black", "grey50"), labels = c("Male", "Female")) +
#	scale_shape_manual(values = c(15, 19), labels = c("Male", "Female"))
#deathPenaltyPlot3

deathPenaltyPlot4 <- ggplot(deathPenaltyPredsDF, aes(y = emmean, x = treatmentGroup, ymin = lower.CL, ymax = upper.CL, color = female, fill = female)) +
	geom_bar(position = position_dodge(), stat = "identity") +
	facet_wrap(~partisanship) + 
	geom_errorbar(position = position_dodge(0.9), width=0.2, color = "grey25") + 
	theme_bw() +
	xlab("\nOpinion Writer Treatment Groups") +
	ylab("Predicted Value - Decision Thermometer\n") +
	ggtitle("Death Penalty Vignette Expanded Results") +
	theme(plot.title = element_text(hjust = 0.5)) +
	scale_x_discrete(labels = c("Control", "Liberal\nMale\nJustice", "Conservative\nMale\nJustice", "Liberal\nFemale\nJustice", "Conservative\nFemale\nJustice")) + 
	scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 10)) +
	theme(legend.position = "bottom", legend.title = element_blank()) +
	scale_color_manual(values = c("grey75", "grey50"), labels = c("Male", "Female")) +
	scale_fill_manual(values = c("grey75", "grey50"), labels = c("Male", "Female"))
deathPenaltyPlot4

#################
### stargazer ###
#################

stargazer(deathPenaltyTreatment, deathPenaltyGender, deathPenaltyPartisanship, modelDeathPenalty, title = "Ordinary Least Squares Regression Results, Decision Thermometer, Death Penalty Vignettes", align = TRUE)

##############################
### CONTROL FOR LEGITIMACY ###
##############################

modelDeathPenaltyLegit <- lm(decisionTherm ~ treatmentGroup
							+ legitScore, 
						data = dataDeathPenalty)
summary(modelDeathPenaltyLegit)
nobs(modelDeathPenaltyLegit)

modelDeathPenaltyLegit2 <- lm(decisionTherm ~ treatmentGroup
							+ legitAverage, 
						data = dataDeathPenalty)
summary(modelDeathPenaltyLegit2)
nobs(modelDeathPenaltyLegit2)

# legit score is a significant predictor

##############################
### BIG MODEL FOR APPENDIX ###
##############################

deathPenaltyWithControls <- lm(decisionTherm ~ treatmentGroup
						+ female
						+ partisanship
						+ income
						+ education
						+ age
						+ female * partisanship
						+ treatmentGroup * female * partisanship
						+ treatmentGroup * female
						+ treatmentGroup * partisanship, 
						data = dataDeathPenalty)
summary(deathPenaltyWithControls)
nobs(deathPenaltyWithControls)

deathPenaltyBigPreds <- emmeans::emmeans(deathPenaltyWithControls, specs = c("female", "partisanship"), by = "treatmentGroup")
pairs(deathPenaltyBigPreds)
deathPenaltyBigPredsDF <- as.data.frame(deathPenaltyBigPreds)

deathPenaltyBigPredsDF$treatmentGroup <- factor(deathPenaltyBigPredsDF$treatmentGroup, levels = c(10, 6, 7, 8, 9))

names <- as_labeller(c(`10`="Control", `6`="Liberal Male Justice", `7`="Conservative Male Justice", `8`="Liberal Female Justice", `9`="Conservative Female Justice"))

deathPenaltyBigPlot <- ggplot(deathPenaltyBigPredsDF, aes(y = emmean, x = partisanship, color = female, shape = female)) +
	geom_point(position = position_dodge(0.3)) + 
	geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0, position = position_dodge(0.3)) +
	facet_grid(~treatmentGroup, labeller = names) + 
	theme_bw() +
	xlab("\nParticipant Partisanship") +
	ylab("Predicted Value - Decision Thermometer\n") +
	ggtitle("Death Penalty Vignette Full Results") +
	theme(plot.title = element_text(hjust = 0.5)) +
	scale_x_discrete(labels = c("Dem", "Ind", "Rep")) + 
	scale_y_continuous(limits = c(0, 125), breaks = seq(0, 125, 10)) +
	theme(legend.position = "bottom", legend.title = element_blank()) +
	scale_color_manual(values = c("black", "grey50"), labels = c("Male", "Female")) +
	scale_shape_manual(values = c(15, 19), labels = c("Male", "Female"))
deathPenaltyBigPlot

stargazer(deathPenaltyTreatment, deathPenaltyGender, deathPenaltyPartisanship, modelDeathPenalty, deathPenaltyWithControls, title = "Ordinary Least Squares Regression Results, Decision Thermometer, Death Penalty Decision", align = TRUE)

############################
### DEATH PENALTY BACKUP ###
############################

dataDeathPenalty$agreeDV <- ifelse(dataDeathPenalty$agreeWithDecisionRecode == 1, 1, 0)
dataDeathPenalty <- dataDeathPenalty %>% filter(!is.na(agreeDV))

modelAgreeDeathPenalty <- glm(agreeDV ~ treatmentGroup
						+ female
						+ partisanship
						+ female * partisanship
						+ treatmentGroup * female * partisanship
						+ treatmentGroup * female
						+ treatmentGroup * partisanship,
						family = binomial(link = logit), 
						data = dataDeathPenalty)
summary(modelAgreeDeathPenalty)
nobs(modelAgreeDeathPenalty)

dpPred2 <- ggpredict(modelAgreeDeathPenalty, terms = c("treatmentGroup [6, 7, 8, 9, 10]", "female [0,1]", "partisanship [Democrat, Republican]", "deathPenaltySupport [1]"))

deathPenaltyPlot2 <- ggplot(dpPred2, aes(y = predicted, x = x, color = group)) +
	geom_point(position = position_dodge(0.3)) + 
	geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, position = position_dodge(0.3)) +
	facet_wrap(~facet) + 
	theme_bw() +
	xlab("\nOpinion Writer Treatment Group") +
	ylab("Predicted Value - Decision Thermometer\n") +
	ggtitle("Death Penalty Vignette") +
	theme(plot.title = element_text(hjust = 0.5)) +
	scale_x_discrete(labels = c("Liberal\nMale", "Conservative\nMale", "Liberal\nFemale", "Conservative\nFemale", "Control")) + 
	scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, 0.1)) +
	theme(legend.position = "bottom", legend.title = element_blank()) +
	scale_color_manual(values = c("black", "azure4"), labels = c("Male", "Female"))
deathPenaltyPlot2

# scTherm as DV
modelDeathPenaltySC <- lm(scTherm ~ treatmentGroup
						+ female
						+ partisanship
						+ female * partisanship
						+ treatmentGroup * female * partisanship
						+ treatmentGroup * female
						+ treatmentGroup * partisanship,
						data = dataDeathPenalty)
summary(modelDeathPenaltySC)
nobs(modelDeathPenaltySC)

dpPred3 <- ggpredict(modelDeathPenaltySC, terms = c("treatmentGroup [6, 7, 8, 9, 10]", "female [0,1]", "partisanship [Democrat, Republican]"))

deathPenaltyPlot3 <- ggplot(dpPred3, aes(y = predicted, x = x, color = group)) +
	geom_point(position = position_dodge(0.3)) + 
	geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, position = position_dodge(0.3)) +
	facet_wrap(~facet) + 
	theme_bw() +
	xlab("\nOpinion Writer Treatment Group") +
	ylab("Predicted Value - Decision Thermometer\n") +
	ggtitle("Death Penalty Vignette") +
	theme(plot.title = element_text(hjust = 0.5)) +
	scale_x_discrete(labels = c("Liberal\nMale", "Conservative\nMale", "Liberal\nFemale", "Conservative\nFemale", "Control")) + 
	scale_y_continuous(limits = c(20, 115), breaks = seq(20, 115, 10)) +
	theme(legend.position = "bottom", legend.title = element_blank()) +
	scale_color_manual(values = c("black", "azure4"), labels = c("Male", "Female"))
deathPenaltyPlot3

